Electron waves in chemically substituted graphene 
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We present exact analytical and numerical results for the electronic spectra and the Friedel oscillations around 
a substitutional impurity atom in a graphene lattice. A chemical dopant in graphene introduces changes in the 
on-site potential as well as in the hopping amplitude. We employ a T-matrix formalism and find that disorder 
in the hopping introduces additional interference terms around the impurity that can be understood in terms of 
bound, semi-bound, and unbound processes for the Dirac electrons. These interference effects can be detected 
by scanning tunneling microscopy. 

PACS numbers: 73.20.Hb,81.05.Uw,73.20.-r, 73.23.-b 
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Introduction. Graphene sl one atom- thick layer of 
graphite, has been intensively studied l'^ S B B F 
II0IIIIIII2I] due to its fascinating physical properties | 
Graphene is a zero-gap semiconductor and its low-energy 
electronic excitations are described in terms of a Dirac spec- 
trum. Because of this, disorder in the form of impurities 
liOH E E4], defects [ll H [17, 18] and surfaces 



i|2d] can have a strong effect, especially when the chemi- 
cal potential crosses the Dirac point. In this work, we consider 
dilute chemical dopants in graphene incorporated as substi- 
tutional atoms, and calculate local single particle properties 
such as the electronic energy spectra, local density of states, 
and electron density distribution, which shows Friedel oscil- 
lations as depicted in Fig. ([T]). The results presented here 
can be measured by scanning tunneling microscopy (STM) 



Atomic substitution in a carbon (C) honeycomb lattice is 
chemically possible for boron (B) and nitrogen (N) atoms. 
There have been several experimental studies of B and N 
substitution in highly-oriented pyrolytic graphite (HOPG) 
graphitic structure s [1261. nanoribbons [27], carbon 
nanotubes ^ 28, 29, 30, lUH 33, 34, 35, 36, 37] and 
fullerenes 113 8n . When B or/and N atoms replace a carbon 
atom in a graphene sheet they have the following effects: (1) 
they act as impurity scattering centers; (2) they act as hole- 
or electron-dopants; (3) they introduce lattice distortions. In 
the case of B there is an increase of the absolute value of the 
hopping amplitude given its larger atomic radius (r 0.85 A) 
as compared to C (r :^ 0.7 A). On the other hand, a N atom 
impurity (r :^ 0.65 A) can be modeled by a smaller hopping 
integral. For the change in the on-site energy, we assume the 
C on-site energy to be zero (our reference state); the smaller 
atomic number of B leads to a positive local energy, whereas 
on a site occupied by N, the local energy should be smaller 
and therefore is modeled as a negative local energy. The ef- 
fect of a different on-site energy has been extensively studied 
in the past [8, 11, 12, 13, 14], but the effect on the hopping has 
been overlooked. In this work we study the combined effect 
of local changes in the atomic energy as well as the hopping 



amplitudes. A vacancy has been modeled by an infinite, on- 
site energy potential, but it is also exactly represented by the 
remotion of particular hopping processes from the Hamilto- 
nian. The opposite limit of very large hopping to the impurity 
site corresponds to a four-atom molecule. 
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FIG. 1: (color on line) Real-space electronic distribution and Friedel 
oscillations around an impurity with ionic radius smaller than carbon 
(such as N) (to = 0.5t). 

Hopping and potential disorder. The honeycomb lattice has 
a unit cell represented in Fig. Q by the vectors ai and a2, 
such that |ai| = \a2\ = a (a = >/3ao — 2.461 A, where ao 
is the carbon-carbon distance). In this basis any lattice vec- 
tor r is represented as r = nai + ma2, with n, m integers. 
In Cartesian coordinates, ai = ao(3, V^, 0)/2 and a2 = 
ao(3, — a/3, 0)/2. The reciprocal lattice vectors are given by: 
61 = 27r(l, V3,0)/(3ao) and 62 = 27r(l, - ^3, 0)/(3ao), 
and the vectors connecting any A atom to its nearest neigh- 
bors are: = (ai — 2a2)/3, ^2 = («2 — 2ai)/3, and 
^3 = {di + ci2)/3. Using these definitions the Hamilto- 
nian can be written as: H = Hq -\- Vt -\- Vo, where Hq = 
-t Er +^'^(^ - a2)a{r)^b\r - ai)a(r) +h.c.], 
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FIG. 2: (color on line) The honeycomb lattice with a carbon atom 
replaced by an impurity atom, such a B or N. The local hopping 
parameter changes from t tot — to. 



is the kinetic energy operator and (b^) are creation opera- 
tors in the A (B) sites. The operators Vt = to[b^{0)a{0) + 
6t(-a2)a(0) + 6t(_ai)a(0) + h.c] and Vo = £oa'^{0)a{0) 
are the impurity terms for hopping and potential disorder, re- 
spectively. In the particular case to = t, the scattering term 
Vt represents a vacancy. In this work we consider the case of 
zero chemical potential when the Fermi level crosses the Dirac 
point. This is the case in which the system is most susceptible 
to the presence of impurities. The equations of motion for the 
Green's functions can be readily written, and are given by: 



iuJnGaa{^n,k,p)=Sk,p-^ 



-2. 
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iuJnG allien, k,p) = 

iu;nGbb{^n,k,p) = Sk,p 
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where Xk,p = -t^p{Sk,p - h/Nct), (^p = 1 + e-'P ''^ + 
g-2p a2^ and Nc is the total number of unit cells in the lattice. 
The sublattice symmetry is broken by the presence of the im- 
purity and this is manifested in the fact that Xpq 7^ Xqp and in 
the asymmetric way in which the eo-tcrm appears in the equa- 
tions above. The above set of equations can be solved exactly. 
The fact that the scattering term Vt depends on (j)^ leads to a 
more complex form for the T-matrix than usual. The exact 
solution for the Green's functions can be written in the form: 



Gaa{k,p) = Sk,pGl + g + h[Gl + Gl]+ GItgI 
G,,(fe,p) = 5u,^Gl + ^-^GlTGf-^ . (2) 



where all the terms also depend on cjn (omitted here for 
brevity). They are defined as: 



giun) = tlG''{u;n)/[NcD{un)], 

hicOn) = to(t-to)/[N,D{uJn)l 



(3) 
(4) 



T{uJn) 



iuJnto{2t - to) - Spt'^ 



(5) 



where the denominator D{(jJn) is given by 



D{iJn) = {t-tof^[iuOnto{2t-to)-eot^]G^{uOn) (6) 

and G^K) = Efc G\ujn, k)/Nc with G^ = G^{cOn. k) = 
{iu;n)/[{i^n)'^ — ^^l^fcp], the diagonal Green's function for 
the clean system. 

The significance of the term g{u;n) in ([B which only 
appears in Gaa, is easily appreciated if we do the double 
Fourier transform to real space. This term only contributes to 
Gaa (0,0), the return amplitude to the impurity site. The factor 
1/D{ujn) contains a sum over an infinite series of intermedi- 
ate scattering events, but the overall process is bounded and 
the tl factor denotes hopping from the impurity to the nearest 
neighbor B- sites and back to the impurity site. A similar in- 
terpretation can be given to another term which only appears 
in Gaa^ namely, G^{u;nj k)h{<jJn)'- this term contributes only 
to Gaai'^j 0) and describes an additional amplitude of propa- 
gation between the impurity site and another A site. No such 
terms can, of course, appear in G55 when the inpurity is at an 
A site. 

The exact general solution for both diagonal components of 
the Green's functions for a substitutional impurity with poten- 
tial and hopping disorder is contained in ([T]) and (EJ) . The lo- 
cal density of states (LDOS) can be obtained from the Green's 
functions (1112b . after analytical continuation icOn co -\- zO+. 
For sites in the A and B sublattices, it is given by: Piy(r, co) = 
— ImG^^(cj, r, r)/7r , where u = a^b and r is the position of 
the unit cell. The local number of electrons, for a half-filled 
band, is obtained by: 



na,b{r) = / dujpa,b{r,uj) 

J-W 



(7) 



where W = 3tis half of the total electronic bandwidth. 

Strong suppression of hopping and vacancies. A vacancy 
corresponds to to = t, leading to a simplified form for Gaa' 



Gaa{k,p) = GlSk,p 
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where Gaa, G^ and G^ also depend on ujn, and the on-site 
impurity potential was set to zero (eo = 0). Because the va- 
cancy creates a three-site zig-zag edge, one expects the ap- 
pearance of zero energy modes [SI]. The l/iujn term in ([5]) 
leads to a contribution, Pa(0, uj) oc 6{uj) in the total DOS, as 
shown in Fig. ([3]). This contribution comes from the atom 
that has been disconnected from the rest of the lattice. It 
corresponds to exactly one state and to properly represent a 
missing atom this contribution to the DOS should be ignored. 
At small energies (uo <C t) we obtain for the nearest neigh- 
bor B-sites that ^55(0, 0, cj) ^ -l/(9t^G°(cc;)), in agreement 



3 



with previous studies ^ [s^ H^] for vacancy modeled as in- 
finite on-site potential, Sq = oo. In this limit, pi){0,uj) o:^ 
(V3|cj|/6^t2){l + (36^727)[cjln|y3cjV6^^'|]"'}givmga 
resonance at = 0. The full numerical calculation for all en- 
ergies is shown in Fig. ^ for the impurity site and in Fig. ^ 
for the nearest neighbor sites, with the spatial dependence of 
the amplitude of the low-energy states shown in the insets. 



atomic number of N atoms gives a negative on-site potential 
£o with respect to the carbon sites. 
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FIG. 3: (Color online) LDOS for to = t and to = 0.9t at the impu- 
rity site. Main graph corresponds to (0, a; ) for to = t (straight blue 
line) and for to = 0.9t (dashed red line). Left inset: pa{r,io = 0) 
for to = t. Right inset: pa{r,uj = -0.075t) for to = 0.9t. Notice 
the scale difference. 
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FIG. 4: (Color online) LDOS for to = t and to = 0.9t at the 
nearest neighbor sites of the impurity. Main graph corresponds to 
Pb{0,co) for to = t (straight blue line) and for to = 0.9t (dashed 
red line). Left inset: pb(r,uj = 0) for to = t. Right inset: 
pt{r,ou = -0.075t) for to = 0.9t. 

For to slightly smaller than t, hopping to the impurity site 
is strongly suppressed but not completely absent. Hybridiza- 
tion between the low energy states at the impurity site and 
at the nearest neighbors sites leads to splitting into two reso- 
nant states, one with energy shifted to a negative value and the 
other to a positive value, as shown in Figs. ^H^. 

Nitrogen substitution. As discussed previously, a N atom 
has a smaller atomic radius than a carbon atom, and this cor- 
responds to a smaller hopping amplitude between the N impu- 
rity and the neighboring carbon atoms. In addition, the larger 



FIG. 5: (color on line) Left: LDOS at the impurity site (above) and 
at its nearest neighbors (below), for to = 0.5t. Right: intensity plots 
of pa{r,u = -0.525t) (top), andpb(r,cj = -0.975t) (bottom). 

We first study the effects of hopping disorder alone, and 
consider the case to = 0.5t. Fig. ([T]) shows an intensity plot 
of the real- space distribution of the number of electrons, given 
by d?]). Interference effects give rise to Friedel oscillations 
displaying the underlying six-fold symmetry. The main con- 
tributions in the sum over negative states ^ come from: (i) 
the resonance states created by the impurity; and (ii) the spec- 
tral weight under the van Hove singularity that is also affected 
by the impurity. Fig. ([5]) shows these two contributions sepa- 
rately. In the left, we show the LDOS spectrum at the impu- 
rity site, and at its nearest-neighbors (B-sites). The spectrum 
shows a resonance peak at ujq = ±0.525t on the impurity, 
while for its nearest-neighbors it is mostly dominated by the 
van Hove peaks at cc; ±t (c^o = ±0.975t). The size of these 
peaks decreases away from the impurity and Fig. (|5]) shows 
the intensity plots of these peaks in real-space, which can be 
directly measured using STM spectroscopy. 

When a negative So is also added, the LDOS gets modi- 
fied as shown in Fig. ^ for the impurity and nearest neigh- 
bor sites. Any finite potential value at the impurity breaks 
particle-hole symmetry and as a result the spectrum becomes 
asymmetric under reflection (uj —uj). 

Boron substitution. In this case the impurity site has larger 
hopping amplitude and a positive on-site potential. Fig. d?]) 
shows the LDOS at the impurity and the nearest neighbor sites 
for a case with to = — 0.5t (so that the hopping amplitude to 
the impurity, t — to, is larger than the homogeneous system) 
and So = 0.525t. The peaks of the LDOS at certain ener- 
gies originate from impurity resonance states and van Hove 
singularities. The real- space intensity plots are at resonance: 
LJo = -2.935tforthe sublattice-A LDOS, andcjo = -0.975t 
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FIG. 6: LDOS for the case to = 0.5t and £o = -0.525t, at the 
impurity site (left) and at the nearest neighbor sites (right). 



atoms have very specific spectroscopic signatures. From the 
knowledge of these Green's function one can also calculate 
the transport properties of chemically substituted graphene. 
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FIG. 7: (color on line) Electronic spectra for the case to = — 0.5t 
and eo = 0.525t, corresponding to an impurity such as boron. Left: 
LDOS at the impurity site (top) and nearest neighbors (bottom). 
Right: real-space intensity plot of pa{r^uj — — 2.935t) (top), and 
pi,{r,uj = -0.975t) (bottom). 

(at the van Hove singularity) for the sublattice-B LDOS. Just 
as for the case of N shown in Fig. ([5]) the van Hove singu- 
larities are more strongly affected at certain directions. The 
van Hove peaks are strongest on B-sites where the LDOS 
can be understood as being originated by three impurities (the 
three nearest neighbor B-sites of the actual impurity site) and 
hence the star-shaped symmetry shown in the intensity maps 
at uoq = — 0.975t in Fig. ([5]) and Fig. (|7]). Near the band edge, 
the band structure of the clean system behave more like a con- 
ventional 2D system, so the resonance peak at uoq = — 2.935t 
has a symmetric spatial distribution (Fig. (|7]), top). 

Conclusions. We find the exact electronic Green's func- 
tions in the presence of an impurity that modifies both the 
local atomic energy as well as the hopping amplitude with 
neighboring C atoms. We have shown that the presence of 
the impurity leads to strong modifications of the local den- 
sity of states and also to the presence of unusual Friedel os- 
cillations, both quantities being accessible to measurement 
through STM spectroscopy. We have applied our theory to 
the case of substitutional B and N and found that these two 
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